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We introduce a new non-Hermitian approximation of Bloch equations. This approximation pro- 
vides a complete description of the excitation, relaxation and decoherence dynamics of single as well 
as ensembles of coupled quantum emitters (atoms or molecules) in weak laser fields, taking into 
account collective effects and dephasing. This is demonstrated by computing the numerical wave 
packet solution of a time-dependent non-Hermitian Maxwell-Schrodinger equation describing the 
interaction of electromagnetic radiation with a quantum (atomic or molecular) nano-structure and 
by comparing the calculated transmission, reflection, and absorption spectra with those obtained 
from the numerical solution of the Maxwell-Liouville-von Neumann equation. We provide the key 
ingredients for easy-to-use implementation of the proposed scheme using time-dependent decay rates 
and identify the limits and error scaling of this approximation. 
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I. INTRODUCTION 

Optics of nanoscale materials has attracted consider- 
able attention in the past several years [1-4] due to vari- 
ous important applications [5]. Exploring electrodynam- 
ics of near-fields associated with subwavelength systems, 
researchers are now truly dwelling into nanoscale [5-9]. 
Owing to both new materials processing techniques [10] 
and continuous progress in laser physics [11] the research 
in nano-optics is currently transitioning from linear sys- 
tems, where materials and their relative arrangement 
control optical properties [12], to the nonlinear regime 
[13]. The latter expands optical control capabilities far 
beyond conventional linear optics as in the case of active 
plasmonic materials [14], for instance, combining highly 
localized electromagnetic (EM) radiation driven by sur- 
face plasmon-polaritons (SPP) with non-linear materi- 
als [15]. 

Yet another promising research direction, namely 
optics of highly coupled exciton-polariton systems, is 
emerging [16, 17]. It basically reincarnates a part 
of research in semiconductors [18, 19], bringing it to 
nanoscale via deposition of ensembles of quantum emit- 
ters (molecules [20-23], quantum dots [24-27]) directly 
on to plasmonic materials. Even in the linear regime, 
when the external EM radiation is not significantly ex- 
citing the quantum sub-system, SPP near-fields can be 
strongly coupled to quantum emitters. This manifests 
itself as a Rabi splitting widely observed in transmission 
experiments [28]. Moreover, a new phenomenon, namely 
collective molecular modes driven by SPP near-fields, has 
been observed [21] and recently explained [29]. It was also 
shown that nanoscale clusters comprised of optically cou- 
pled quantum emitters exhibit collective scattering and 
absorption [30] similar to Dicke superradiance [31]. It 
is hence important to be able to account for collective 
effects in a self-consistent manner. We also note that a 



series of works by Neuhauser et al. [32-34] clearly demon- 
strated that even a single molecule in a close proximity 
of a plasmonic material can significantly alter scattering 
spectra. 

In many applications (such as optics of molecular lay- 
ers coupled to plasmonic materials [20-23, 29], for in- 
stance) self-consistent modeling relies heavily on the nu- 
merical integration of the corresponding Maxwell-Bloch 
equations assuming that static emitter-emitter interac- 
tions can be neglected (which is true for systems at rel- 
atively low densities). Such an approximation results in 
expressing the local polarization in terms of a product of 
the local density of quantum emitters and the local aver- 
aged single emitter's dipole moment [35] . One of the first 
efficient numerical schemes for simulations of nonlinear 
optical phenomena of quantum media driven by exter- 
nal classical EM radiation was proposed by Ziolkowski 
et al. [36]. Using a one-dimensional example of ensem- 
bles of two-level atoms it was shown that the correspond- 
ing Maxwell-Bloch equations can be successfully inte- 
grated using an iterative scheme based on the predictor- 
corrector method [strongly coupled method). Later on 
this approach has been extended to two- [37] and three- 
dimensional systems [38] . Although such a scheme accu- 
rately captures the system's dynamics, it can become ex- 
tensively slow for multidimensional systems [30]. More- 
over this method is limited to two-level systems only. A 
more efficient technique based on the decoupling of Bloch 
equations from the Ampere law was proposed in 2003 by 
Bidegaray [39]. This latter method, usually referred to 
as a weakly coupled method, noticeably improves the ef- 
ficiency of the numerical integration of Maxwell-Bloch 
equations, allowing to consider multilevel quantum me- 
dia [30]. 

The approach proposed in the present paper is based 
on this weakly coupled method [30]. It further improves 
the numerical efficiency of the Maxwell-Bloch integrator 
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for ensembles of multilevel quantum emitters. By incor- 
porating a new non-Hermitian wave packet propagation 
technique into the weakly coupled method, we demon- 
strate that our approach can be successfully applied to 
ensembles of multilevel atoms and diatomic molecules. 

The paper is organized as follows. We first introduce 
our model in section II A, using, as an example, an ensem- 
ble of two-level atoms. The non-Hermitian wave packet 
approximation is then described in section II B. Applica- 
tions of the proposed method to the case of a nano-layer 
comprised of two-level atoms are discussed in section II C. 
We then generalize our method to the case of interact- 
ing multi-level emitters in section II D. Finally, using this 
generalized approach, we calculate in section II E the op- 
tical properties of a nano-layer of coupled molecules, tak- 
ing into account both the vibrational and rotational de- 
grees of freedom, and we reveal several interesting fea- 
tures in the absorption spectra. Last section III summa- 
rizes our work. 



II. THEORETICAL MODEL AND 
APPLICATIONS 



As pointed out in the previous section, the importance 
of collective effects manifested in recent experiments 
and theoretical papers calls for the development of self- 
consistent models capable of taking into account mutual 
EM interaction in ensembles of quantum emitters. While 
direct numerical integration of corresponding Maxwell- 
Bloch equations can be based on either a strong coupled 
method [36] or a weakly coupled one [39] , for multi-level 
systems in two- and especially three-dimensions such a 
brute-force approach becomes numerically very expen- 
sive, both in terms of computation (CPU) times and 
memory requirements. Indeed, the main disadvantage of 
such approaches is that both the CPU times and mem- 
ory requirements scale generally at least as N?, where 
Nh is the dimension of the Hilbert space of the system. 
This unfavorable scaling is directly related to the size of 
the reduced density matrix used to describe the system 
and to the associated number of time-dependent equa- 
tions one has to solve to follow the system's dynamics. 
By contrast, one has only to solve a reduced set of Nh 
time-dependent equations when the system can be de- 
scribed by a single wave function. It would therefore 
be extremely useful to be able to derive a "Schrodinger- 
type" approximation which could describe on an equal 
footing the field-induced coherent dynamics of a multi- 
level quantum system and the associated relaxation and 
decoherence processes. This goal has been achieved in 
the past using different approaches, the three main con- 
tributions being the stochastic Schrodinger method used 
in conjunction with a Monte Carlo integrator [40-45], the 
Gadzuk jumping wave packet scheme [46, 47] and the 
variational wave packet method [48, 49]. These methods 
require the propagation of Nf different wave functions 
and are therefore mainly attractive when Nf <C Nh- 



The approach we propose here is different since it is 
based on perturbation theory. It also significantly speeds 
up calculations since it only requires the propagation of 
a single wave function under the action of an easy-to- 
implcment time-dependent effective Hamiltonian in order 
to reproduce accurately the full dynamics. We demon- 
strate the efficiency and accuracy of our method using, 
first, ensembles of two-level atoms, and then, interacting 
diatomic molecules including both the vibrational and 
rotational degrees of freedom. We also discuss the con- 
ditions under which our method is no longer valid. All 
results are compared with data obtained via direct inte- 
gration of Maxwell-Bloch equations using a weakly cou- 
pled method [30]. 



A. Atomic two-level system 

Let us first consider a two-level quantum system in- 
teracting with EM radiation. We label the two levels as 
|0) and |1), with associated energy eigenvalues Huq and 
huj\, respectively. The corresponding density matrix p(t) 
satisfies the dissipative Liouville-von Neumann equation 
[501 
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where H — Hq + V{ (t) is the total Hamiltonian and f is a 
superoperator, taken in the Lindblad form [51], describ- 
ing relaxation and dephasing processes under Markov ap- 
proximation. The field free Hamiltonian reads 

ff = ^ |0)(0|+^ 1 |l)(l|. (2) 

The interaction of the two-level system with EM radia- 
tion is taken in the form 



Vi{t) = Ml{t){\\) 



(II 



(3) 



£l(t) denotes here the instantaneous Rabi frequency as- 
sociated with the coupling between the quantum system 
and an external field. In Eq. (1), the non-diagonal ele- 
ments of the operator T include a pure dephasing rate 
7* , and the diagonal elements of this operator consist of 
the radiationless decay rate V of the excited state. Equa- 
tions (l)-(3) lead to the well-known Bloch optical equa- 
tions [35] describing the quantum dynamics of a coupled 
two-level system 
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where Wg = u>i — co is the Bohr transition frequency and 
where the dot denotes the time derivative. 
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We assume in the following that the system is initially 
in the ground state |0), and we will show that, under some 
assumptions, the subsequent induced excitation and re- 
laxation dynamics can be accurately described by a non- 
Hermitian wave packet approximation. 



B. Non-Hermitian two-level wave packet 
approximation 

Within the aforementioned simplified model, the sys- 
tem's wave packet |^(t)) can be expanded as 



sity matrix p s (t) 



|*(t)) = co(t)|0)+ Cl (i)|l> 



(5) 



Let us now assume that the ground and excited states 
energies include an imaginary part that we denote as 
+fryo/2 and —h'y 1 /2, respectively. Inserting expan- 
sion (5) into the time-dependent Schrodinger equation 



(6) 



and projecting it onto states |0) and |1) yields the fol- 
lowing set of coupled equations for the coefficients c„(i) 
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We will now derive the differential equations describ- 
ing the temporal dynamics of the products pfj = Cj c*j , 
where the subscript s corresponds to this simplified non- 
Hermitian Schrddinger-type model. Our goal is to ob- 
tain the decay rates 70 and 71 which will allow for an 
approximate description of the system's excitation and 
relaxation dynamics. For the evolution of the popula- 
tions Po (t) and pli(t), one gets 



Poo = (Poi - P10) + 7o Poo 

Pil = ^(*) (pio - P01) - 7i P11 • 

The conservation of the total norm 



thus results in 



Poo(t)+Pu(t) = 



7oPoo(*) = 7i Pu(t) 



(8a) 
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It therefore appears that, since the populations Poo(i) 
and pli(t) are generally time dependent, at least one of 
the two rates 70 or 71, which are not yet fully defined, 
must be taken as a time- dependent function. 

Taking Eq. (10) into account, one finally obtains the 
following set of Bloch equations for the approximate den- 
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to be compared with the exact equations (4a)-(4d). 

Two obvious choices can then be made for the empiri- 
cal decay parameters 70 and 71: 

• The first choice Ji = T allows to reproduce cor- 
rectly the equations (4a) and (4d) describing the 
populations at the cost of degrading the descrip- 
tion of the coherences poi(^) an d Pio{t)- 

• The second choice, 71—70 = 27* + V, allows to 
reproduce correctly the equations (4b) and (4c) de- 
scribing the coherences, at the cost of degrading 
the description of the populations. 

The optimal choice for applications in linear optics of 
nano-materials should be as follows: in weak fields, the 
variation of the populations, as described by perturba- 
tion theory, is a second order term with respect to the 
coupling amplitude while the variation of the coherences 
is a first order term. For a correct description of the 
quantum dynamics in weak fields, it is important to de- 
scribe first order terms accurately. Therefore we proceed 
with the second choice 



71 - 70 = 27* + r. 



(12) 



From Eqs.(lO) and (12) we obtain a set of two empirical 
time dependent decay rates 70 (t) and 71 (t) that we can 
insert into the Schrodinger-type approximation 
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We finally obtain the following set of time-dependent cou- 
pled equations 
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This system is solved numerically using a fourth order 
Runge-Kutta algorithm [52]. Eqs. (14a), (14b) include 
two non-linear non-Hermitian terms which can (as we will 
illustrate below) accurately reproduce the dissipative dy- 
namics in weak fields, i.e. for |ci| 2 <?C | Co | 2 ~ 1- Indeed, 
as one can easily show that, in this limit, Eqs. (lla)-(lld) 
are strictly equivalent to the exact Bloch Eqs. (4a)-(4d). 
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C. Application to a uniform nano-layer of atoms 

While our approximation is equally applicable to three- 
dimensional systems, we have chosen, for the sake of 
simplicity, to test it on a simplified one-dimensional sys- 
tem consisting of a uniform infinite layer of atoms whose 
thickness Az lies in the range of a few hundred nanome- 
ters. An incident radiation field propagating in the pos- 
itive z-direction is represented by a transverse-electric 
mode with respect to the propagation axis. It is char- 
acterized by a single in-plane electric field component 
E x (z,t) and a single in-plane magnetic field component 
H y (z, t). To account for the symmetry of the atomic po- 
larization response, the atoms in the layer are described 
as two-level systems with the following states: an s-type 
ground state and a p x -type excited state. This model 
is a one-dimensional simplification of a more general ap- 
proach used in Ref. [30] and we refer the reader to the 
body of this paper for details. 

The time-domain Maxwell's equations for the dynam- 
ics of the electromagnetic fields 

Hod t Hy = -d z E x (15a) 
e d t E x = -d z H v -d t P x (15b) 

are solved using a generalized finite-difference time- 
domain technique where both the electric and magnetic 
fields are propagated in discretized time and space [53]. 
In these equations, po and eo denote the magnetic per- 
meability and dielectric permittivity of free space. The 
macroscopic polarization of the atomic system 

P x (z,t) = n(fi x ) (16) 

is taken as the expectation value of the atomic transition 
dipole moment p X) where n is the atomic density. 

A self-consistent model is based on the numerical inte- 
gration of Maxwell's equations (15a) and (15b), coupled 
via Eq. (16) to the quantum dynamics. In the mean-field 
approximation employed here it is assumed that the den- 
sity matrix of the ensemble is expressed as a product of 
density matrices of individual quantum emitters driven 
by a local EM field. In order to account for dipole-dipole 
interactions within a single grid cell we follow [54] and 
introduce Lorentz-Lorenz correction for a local electric 
field term into quantum dynamics according to 

p 

E x local = E x + - — , (17) 

3e 

where E x is the solution of Maxwell's equations (15a), 
(15b) and macroscopic polarization is evaluated accord- 
ing to Eq. (16). 

The quantum dynamics is evaluated by computing 
the atomic dipole moment either using the single- 
atom density matrix (4a)- (4d) or the single-atom wave 
packet (14a)-(14b). To compare two approaches in the 
linear regime (i.e. for p n <C poo ~ 1, |ci| 2 < |c | 2 ~ 1), 
we calculate the transmission T(E), reflection R(E) and 



absorption A(E) spectra of an atomic layer of thickness 
Az = 400 nm as a function of the incident photon energy 
E for an atomic transition energy of Eg = Hlub = 2eV. 
The results are shown in Fig. 1. 

The transmission T(E) and reflection R(E) spectra are 
obtained from the normalized Poynting vector 





E X Hy 






inc 



at a specific location under and above the atomic layer, 
respectivelly, where E x , H y and E Xj i nc , H y ^ nc are the 
Fourier components of the total and incident EM fields. 
The absorption spectrum is then simply obtained as 

A(E) = 1 - T(E) - R{E) . (19) 

In Fig. 1 the solution of Maxwell-Liouville-von Neu- 
mann equations is shown as a blue solid line while the 
open red squares are from the solution of our approxi- 
mate non-Hermitian Schrodinger model. One can notice 
a perfect agreement of the two methods irrespective of 
the atomic density. These results clearly demonstrate 
that, in the linear regime where the atomic excitation 
probability remains small and varies linearly with the in- 
cident field intensity, a simple wave packet propagation is 
sufficient to mimic the excitation and relaxation dynam- 
ics of an ensemble of quantum emitters. Therefore in 
the weak field regime, the propagation of the full density 
matrix is unnecessary. 

At low density, most of the incident radiation is sim- 
ply transmitted and a small part of the incident energy 
is absorbed by the atomic ensemble at energies close 
to the atomic transition energy of 2 eV. The absorp- 
tion spectrum shows the conventional Lorentzian profile. 
At higher densities, the absorption spectrum is strongly 
modified due to the appearance of collective excitation 
modes [30] , and a large part of the incident energy is ei- 
ther absorbed or reflected from the atomic nano-layer at 
photon energies close to the transition energy. 

The blue lines with squares in Fig. 2 show, as a 
function of the incident field intensity, the relative error 
\As(Eb) — A L (E B ) \ I A L (E B ) obtained in the calculation 
of the absorption spectrum at the transition energy Eb 
using the Schrodinger approximation As(Eg) when com- 
pared to the solution of the full Liouville-von Neumann 
equation Al(Eb)- The solid blue line is for the lowest 
atomic density n — 2.5 x 10 25 m~ 3 and the dashed blue 
line is for the highest density n — 2.5 x 10 27 m~ 3 . 

One can see that irrespective of the atomic density 
the relative error of the non-Hermitian wave packet ap- 
proximation scales linearly with the field intensity. This 
is not surprising. Indeed, solving coupled Schrodinger 
equations (14a)- (14b) is strictly equivalent to solving the 
approximate Bloch equations (lla)-(lld). The latter dif- 
fer from the exact Bloch equations (4a)- (4d) by a term 
proportional to the excited state population pn. The 
maximum excited state population /0™ ax is also shown 
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FIG. 1: (Color online) Transmission T(E) (panels a,d,g), re- 
flection R(E) (panels b,e,h) and absorption A(E) (panels c,f,i) 
spectra of an atomic layer of thickness Az — 400 nm as a func- 
tion of the incident photon energy E. The atomic density is 
ti = 2.5 x 10 25 m~ 3 in the first column, n = 2.5 x 10 26 m" 3 in 
the second column, and n = 2.5 x 10 27 m -3 in the last column. 
The decay rate and pure dephasing rate are V = 10 12 s _1 and 
7* = 10 15 s _1 , respectively. The atomic transition energy is 
hws = 2eV and the transition dipole moment is 2D. The so- 
lution of Maxwell-Liouville-von Neumann equations is shown 
as a blue solid line while the open red squares are from the so- 
lution of our approximate non-Hermitian Schrodinger model. 




I (a.u.) 

FIG. 2: (Color online) Log-Log plot of the maximum excited 
state population (red solid line with circles) and the relative 
error (blue lines with squares) in the calculation of the ab- 
sorption spectrum A(E) at the transition energy Eb = Hub 
using the Schrodinger approximation when compared to the 
solution of the full Liouville-von Neumann equation as a func- 
tion of the incident field intensity in atomic units. The blue 
solid line is for the lowest atomic density n — 2.5 x 10 25 m~ 3 
and the blue dashed line is for the highest atomic density 
n — 2.5 x 10 27 m -3 . All other parameters are as in Fig. 1. 



in Fig. 2 as a function of the field intensity. It can be 
seen that it also varies linearly in the present weak field 
regime. We can conclude (from Fig. 2 and numerous 
calculations we have performed) that as long as the ex- 
cited state population remains smaller than 1% our wave 
packet approximation can be used safely. 

D. Generalization to multi-level systems 

In order to generalize our Schrodinger-type approxi- 
mation of the excitation and dissipation dynamics of an 
ensemble of quantum emitters to a multi-level system 
we now consider the case of neutral diatomic molecules. 
More specifically, we have chosen the particular case of 
the ground and first excited electronic states of the Li2 
molecule and we follow the electronic dynamics and the 
nuclear motion by expanding the total molecular wave 
function ^{r e , R, t) using the Born-Oppenheimer expan- 
sion 

V{r e ,R,t) = X g(R,t)$ g (r e \R)+Xe{R,t)$ e (r e \R) (20) 

where $> g {r e \R) and $ e (r e \R) denote the electronic wave 
functions associated with the ground X and first 

excited A electronic states of Li2, respectively. 



The electron coordinates are denoted by the vector r e , 
and the vector R = (R,R) represents the internuclear 
vector. 

We now separate the global electronic coordinate r e 
of all electrons into the coordinate r c of the core elec- 
trons and the coordinate r of the active electron [55]. 
The ground X electronic state is considered as a 

2scr state, and the electronic wave function § g (r e \R) is 
expressed in the molecular frame (Hund's case (b) repre- 
sentation) as the product 

$ g (r e \R) = <Mr c \R) R x (r) Y 00 (f) (21) 

where Rx{r) and loo(^) are the radial and angular parts 
of the electronic wave function associated with the ac- 
tive electron. Similarly, the 2pcr excited state of A 
symmetry is expressed in the molecular frame as 

$ e {r e \R) = ct> e {r c ,r\R)R A {r)Y 10 {r). (22) 

Due to the E symmetry of both electronic states, the 
ro- vibrational time-dependent wave functions \ g (R,t) 
and Xe(Rit) can be expanded on a limited set of nor- 
malized Wigner rotation matrices in order to take into 
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account the rotational degree of freedom, following [56] where c g (t) and c v (t) are time-dependent complex coef- 
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Xe (R,t) = }^X e N,M( R ^ D M,o( R )' ( 23b ) 

N,M 

where N denotes the molecular rotational quantum num- 
ber while M denotes its projection on the electric field 
polarization axis x of the laboratory frame. 

Introducing these expansions in the time-dependent 
Schrodinger equation (6) describing the molecule-field 
interaction and projecting onto the electronic and rota- 
tional basis functions yields, in the dipole approximation, 
the following set of coupled differential equations for the 
nuclear wave packets \% m( R > an< ^ Xn m( R i 
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where pax(R) = (RA\r\Rx)r is the electronic transition 
dipole. The ro- vibrational nuclear Hamiltonians li 9 J e (R) 
are defined as 
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where p denotes the molecular reduced mass. The po- 
tential energy curves V g (R) and V e (R) associated with 
the ground and first excited electronic states of the 
molecule are taken from Ref. [571 and the matrix ele- 



ments M. 
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which couple the nuclear wave packets 



evolving on these electronic potential curves can be writ- 
ten using 3j-symbols as 
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We assume that the molecules are prepared at time 
t = in the ro-vibrational level v = and N = of the 
ground electronic state. In weak linearly polarized fields 
and except for a phase factor this ground state compo- 
nent of the molecular wave function remains unaffected 
and the excited state component is limited to N = 1 
and M = 0. The ground and excited nuclear wave pack- 
ets are thus finally expanded in terms of ro-vibrational 
eigenstates as 



X 9 0fi (R,t) = c g (t)ip 3 Qfi (R) 



(27a) 
(27b) 



ficients. (p 9 J^{R) denote here the bound ro-vibrational 
eigenstates of the g/e electronic potentials [58]. It is 
not necessary here to take into account the dissociative 
nuclear eigenstates associated with the excited potential 
since their coupling with the ground vibrational level of 
the ground electronic state is negligible. 

We thus arrive at a multi-level system which is very 
similar to the atomic case described in sections II A and 
II B, except that the initial ground state is now coupled 
with a large set of excited levels. For convenience and for 
an easy comparison with the atomic case, we will label 
the ground state as state number and the excited states 
as states number j ^ 1. Our reference calculations will 
be based on the numerical solutions of the corresponding 
Liouville-von Neumann equations 
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where Hujj is the total energy of the excited state j and 
where Oj(i) denotes the instantaneous Rabi frequency 
associated with the molecule-field interaction as defined 
in Eqs. (24a) and (24b). 7* and T denote the pure de- 
phasing rate and the relaxation rate associated with the 
excited states, respectively. 

In comparison with this "exact" model, our Schrodin- 
ger-type approximation will be based on the numerical 
solutions of the coupled equations for the time-dependent 
expansion coefficients 



i c 



■ 7o(t) 



CO 



E^W c j 



i Cj = flj (t) Co + 



• 7j(*) 



(29a) 
(29b) 



where the time-dependent decay rates 70 (t) and jj(t) are 
now defined as 



7o (t) 



(27* +r)E 3>1 |cj(*)| 2 
|cd(t)| a -£j*iM*)l a 

(27*+r)|c (Q| 2 
|co(i)| 2 -E,>r|ciW| 2 



(30) 
(31) 



One can show that with such a definition of 70 (t) 
and 7j(i), Eqs. (29a) and (29b) are strictly equivalent to 
Eqs. (28a)-(28d) in the limit of weak couplings. In the 
next section we demonstrate that our Schrodinger-type 
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n = 2.5 10 25 m 3 n = 2.5 10 27 m" 




E (eV) E (eV) 



FIG. 3: (Color online) Absorption spectra A(E) of a L12 
molecular layer of thickness Az = 400 nm as a function of 
the incident photon energy E. The molecular density is n = 
2.5 x 10 25 m -3 in the left column (panels (a) and (b)) and 
n — 2.5 x 10 27 m -3 in the right column (panels (c) and (d)). 
The solutions of Maxwell-Liouville-von Neumann equations 
are shown as blue solid lines in the first raw (panels (a) and 

(c) ) while the red solid lines (inverted spectra, panels (b) and 

(d) ) are from the solutions of our approximate non-Hermitian 
Schrodinger model. All other parameters are as in Fig. 1. 



model accurately reproduces the excitation and dissipa- 
tion dynamics of multi-level quantum systems in the limit 
of weak couplings, i.e. for Yl \ c j\ 2 l c o| 2 ~ L 

E. Application to a uniform nano-layer of 
molecules 

We consider a simplified one-dimensional system sim- 
ilar to the one discussed in case of two-level atoms. A 
uniform infinite layer with a thickness of Az — 400 nm 
comprised of Li2 molecules is exposed to incident lin- 
early polarized radiation. The incident field propagates 
in the positive z-direction. We calculate the absorption 
spectrum A(E) of this molecular layer just as we did in 
section II C for atoms. 

To ascertain the validity of the proposed Schrodinger- 
type approximation, absorption spectra are represented 
in Fig. 3 as a function of the incident photon energy 
E, in the linear regime, for two different molecular den- 
sities: n — 2.5 x 10 25 m~ 3 in the left column (panels 

(a) and (b)) and n = 2.5 x 10 27 m~ 3 in the right column 
(panels (c) and (d)). The solutions obtained via integrat- 
ing Maxwell-Liouville-von Neumann equations are shown 
in the upper panels (a) and (c) as blue solid lines while 
the ones obtained from our approximate non-Hermitian 
Schrodinger model are shown in red in the lower panels 

(b) and (d) (inverted spectra). 



As in the case of two- level atoms, one can notice a 
perfect agreement of the two methods irrespective of the 
molecular density. This shows that, in the linear regime 
where the molecular excitation probability remains small 
and varies linearly with the incident field intensity, a sim- 
ple wave packet propagation is sufficient to mimic the 
excitation, dissipation and decoherence dynamics of an 
ensemble of multi-level quantum emitters. The propaga- 
tion of the full density matrix is then, again, unnecessary. 

It is important to note that at low density, we observe 
a series of overlapping vibrational resonances which re- 
flects the vibrational structure of the excited molecular 
potential and which follows the Franck-Condon principle 
[59, 60]. The green labels seen in panel (a) of Fig. 3 in- 
dicate the excited state vibrational level responsible for 
the observed resonance. 

At higher density, the absorption spectrum is strongly 
distorted and one observes, just like in the atomic case 
[30], the appearance of collective excitation modes, the 
difference being that a vibrational structure may still be 
present in some of these collective molecular excitation 
modes. A detailed analysis of the physics underlying the 
appearance of these intriguing molecular collective modes 
will be presented in another paper. For high densities, 
we could also observe that the radiation field is not trans- 
mitted anymore through the molecular layer in the ab- 
sorption window 1.2 eV ^ E ^ 2.5 eV of the molecule. 
Within this window, the field is either absorbed or re- 
flected. The small inset seen in panel (c) of Fig. 3 shows 
a magnification of the absorption spectrum in the energy 
range 1.7 eV ^ E ^ 2.1 eV corresponding to the "nor- 
mal" low-density aborption spectrum. One can see in 
this inset, and by comparing with panel (a) of Fig. 3, 
that the absorption spectrum is not strongly modified at 
high molecular density in this energy region. 

Figure (4) finally shows the CPU time necessary for the 
calculation of these absorption spectra as a function of 
the number of quantum states introduced in this multi- 
level model for both the Maxwell-Liouville-von Neumann 
(blue line with circles) and Maxwell-Schrodinger (red line 
with squares) approaches. One can observe a substan- 
tial difference in CPU times which can prove of crucial 
importance when one has to deal with realistic three- 
dimensional systems. This origin of the observed gain in 
CPU time relies on the necessity of propagating a single 
wave function instead of a full density matrix. 



III. SUMMARY 

We proposed a new non-Hermitian approximation of 
Bloch optical equations. Our method provides an accu- 
rate, complete description of the excitation, relaxation 
and decoherence dynamics of single as well as ensem- 
bles of coupled quantum emitters (atoms or molecules) 
in weak EM fields, taking into account collective effects 
and dephasing. We demonstrated the applicability of the 
method by computing optical properties of thin layers 
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FIG. 4: (Color online) CPU process time necessary on a 
Intel Xeon E5-1650 processor for the calculation of the ab- 
sorption spectrum of a Li2 molecular nano-layer of thickness 
Az — 400 nm as a function of the number of quantum levels 
included in the calculation. The blue line with circles is for 
the solution obtained from Maxwell-Liouville-von Neumann 
equations while the red line with squares is for our proposed 
Schrodinger-type approximation. The spatial grid has a total 
size of 2560 nm with a spatial step of 1 nm. The time prop- 
agation is performed on a temporal grid of total size 1.7 ps 
with a time step of 1.7 as. 



comprised of two-level atoms and diatomic molecules. It 
was shown that, in the limit of weak incident fields, the 
dynamics of interacting quantum emitters can be suc- 
cessfully described by our set of approximated equations, 
which result in a substantial gain both in CPU time and 
computer memory requirements. These calculations also 
reveal some intriguing new collective molecular excita- 
tion modes which will be presented in detail in another 
publication. The proposed approach was demonstrated 
to provide a substantial increase in numerical efficiency 
for self-consistent simulations. 
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